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Abstract. - We investigate the coalescence of two DNA- bubbles initially located at weak seg- 
ments and separated by a more stable barrier region in a designed construct of double-stranded 
DNA. The characteristic time for bubble coalescence and the corresponding distribution are 
derived, as well as the distribution of coalescence positions along the barrier. Below the melting 
temperature, we find a Kramers-type barrier crossing behaviour, while at high temperatures, 
the bubble corners perform drift-diffusion towards coalescence. The results are obtained by 
mapping the bubble dynamics on the problem of two vicious walkers in opposite potentials. 



Introduction. - The Watson-Crick double helix is the thermodynamically stable state 
of double-stranded DNA in a wide range of temperatures and salt conditions [1]. This sta- 
bility is effected by hydrogen bonds between the bases in individual base-pairs (bps), and 
the stronger stacking interactions between nearest neighbour pairs of bps. Driven by thermal 
fluctuations double-stranded DNA can break apart, to form denaturation bubbles of flexible 
single-stranded DNA [2]. Although rare, bp-opening events expose active groups of DNA 
bases, that are otherwise buried within the double helix. They are crucial for the interaction 
with proteins and chemicals, and therefore for the biological function of DNA. Bubble kinetics 
has been probed in NMR studies [3] , and the growth and reannealing dynamics of individual 
bubbles has been measured in real time in single DNA fluorescence correlation setups [4] . 

© EDP Sciences 



2 



EUROPHYSICS LETTERS 



61) 



I 



2k 



63) 



2k. 



1 1 1 1 1 

I I I I L 



II I hy 

Mil /*" 



I I I I T 

1 1 1 1 L 





I I I I T 



TTT 
J_I_L 



I I I 
1 1 1 



X 

TTT 



' ' ' 



I I I I 
'''' 



^TTT 



Y N 
I I I I 



' ' ' I 



/ T <T.<T b \ 



1 



1 



/ \T, < T\ 



time 




Fig. 1 - Left: Schematic of the bubble coalescence setup in a designed DNA construct. It is clamped 
at both ends and consists of two outer soft zones (thin red lines) of lengths Nl , Nr bps with melting 
temperature T 3 and a stronger iV-bps-long barrier zone (thick blue lines) with TJ, > T s . a) All bps 
closed (T < T 3 < Tb). b) Soft zones open by raising the temperature above T a . bi-bs) Successive 
opening of the barrier driven mainly by fluctuations (T < Tb) or drift (T > Tb) until coalescence. The 
discrete coordinates X, Y = 0, . . . , N are defined as the positions of the interfaces between the closed 
and broken bps. Right: Plot of the linear potentials experienced by the respective bubble interfaces 
in the case T < T b (f < 0) in terms of the dimensionless quantities x, y, f (see text). 



The delicate sensitivity of bubble dynamics and, therefore, bubble nucleation to the local 
DNA bp-sequence [5,6] suggests a new method to use single molecule tools to obtain indepen- 
dently DNA stability parameters, as sketched in figure^ A short stretch of DNA, clamped at 
both ends, is designed such that two soft zones consisting of weaker AT-bps are separated by 
a more stable barrier region rich in GC bps. For simplicity, we assume that both soft zones 
and barrier are homopolymcrs with a bp-dissociation free energy s' and e, respectively, and, 
in accordance with the experimental findings of reference [4] , we neglect secondary structure 
formation in the barrier zone. At temperatures higher than the melting temperature T s of 
the soft zones but still lower than the melting temperature Tb of the barrier region, thermal 
fluctuations will gradually dissociate the barrier, until the two bubbles coalesce. We also 
study the case when the system is prepared as above and then T suddenly increased such that 
T > Tb > T s so that the system is driven towards coalescence. In both cases the two bound- 
aries between bubbles and barrier perform a (biased) random walk in opposite free energy 
potentials. The quantities of interest are the bubble coalescence time and position. It turns 
out that this problem can be mapped on a previously unsolved case of two vicious walkers in 
opposite potentials in one dimension and, therefore, is of interest by itself. 

Formulation of the problem. - The statistical weight of the construct before coalescence 
3f X Y = ^ e N ^'y(x-Y+N)[3e^ e N R f3e^ & t T b > T > T s involves a coopcrativity (ring) factor 
£ w 10~ 5 for each bubble, and a Boltzmann factor for each broken bp with free energies 
e' > and e < 0, compare reference [B]^ 1 ) Upon coalescence, the boundary free energy 
corresponding to one factor £ is released, 3? C oa\ = £e( NL+NR ^ e + Ar / 3e ; stabilizing the system 
against immediate transition back to a two-bubble state. It should therefore be possible to 
experimentally detect the coalescence. In our theoretical approach, this corresponds to an 
absorbing boundary for the random walking interfaces at X = Y. 



( 1 )We denote a stable state by a negative free energy. 
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Fig. 2 - Stochastic trajectories of the random walk of the two fork positions encompassing the double- 
helical bridge of N = 50 bps between the vicinal bubbles for T s < T < Tb with u — exp(/3e) = 0.98 
(left) and T > TJ, > T s , u = 1.10 (right). Note that below the melting temperature TJ, of the barrier 
zone the trajectories have the tendency to move toward the reflecting boundaries at the corners of the 
soft zone and the coalescence takes longer time due to the presence of the potential barrier. On the 
other hand, above Tb the trajectories show a tendency to move fast toward the centre of the barrier 
zone, as expected from the funnelling nature of the positive force directed toward the middle. 



The dynamics of this, intrinsically discrete, system can be modelled by methods developed 
previously, namely, the stochastic Gillespie scheme [7] , which allows one to obtain single bubble 
trajectories and thereby the access to the noise relevant in single molecule experiments, and 
the master equation approach [6,8,9]. As detailed in reference [10], the explicit sequence of 
bps is included into these approaches. Typical examples of individual trajectories resulting 
from the Gillespie scheme are displayed in figure |3 where traces of the two interfaces ("forks") 
cornering the barrier region are shown. Bubble coalescence terminates each pair of trajectories. 

Both the Gillespie scheme and the master equation are based on a specific choice for 
the rate constants, t^~(Q), of (un)zipping the bp at the position Q = A, Y of one of the 
zipping forks where double-stranded DNA branches into two strands of single-stranded DNA. 
As discussed in detail in references [6,8-10], we define 

tr = k/2; t+(Q) = ku(Q)s(m)/2, s(m) = {(m + l)/(m + 2)} c , (1) 

together with the appropriate boundary conditions. Here, t - is the position-independent 
zipping rate k/2 for a bp at the zipping fork, while the bp-unzipping at location Q involves 
the Boltzmann factor u(Q) = cxp((3e(Q)) for disrupting the bp at Q, as well as the factor s(m) 
that stems from the entropy loss of forming a closed polymer ring. The s-factor depends on 
the size m of the bubble cornered by the respective bp, and c m 1.76 is the critical exponent. 
In our illustrations and analytic approach we use the reflecting boundary conditions at the 
edges of the barrier regions, i.e., we assume j3e' 3> 1 so that the soft zones are always open. 

The rates t define the transfer matrix W, that governs the random walk of the zipping 
forks at either end of the barrier region. The coalescence dynamics is then quantified by the 
probability distribution Pjj [X, Y, t) to find the left and right zipping forks at positions X and 
Y, as controlled by the master equation dPo{X,Y,t) / dt = WPd (X, Y, t) . Solution of this 
equation by cither method, as detailed elsewhere [6-8, 10], yields the quantities of interest 
such as the mean coalescence time or distribution of the coalescence position. 
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Continuous description and semi-analytic solution. Both the master equation and 
Gillespie approaches can be used for arbitrary bp sequences and their usage is straightforward, 
however, they both have their limitations. The Gillespie algorithm performs badly in the case 
of a relatively strong barrier (/3e <C —1) when the coalescence time grows exponentially. 
This demands exponentially increasing simulation time. The master equation is becoming 
numerically prohibitive with increasing size of the barrier. On a regular PC the length of 
about N w 100 bps is close to the numerical limit of our master equation approach. 

For the designed construct we can find an efficient semi-analytical solution for its dynamics. 
Due to the large length of the soft zones Nl,r > lwe can well approximate the bubble entropic 
factors s(m) = {(m + l)/(m + 2)} c by 1. Furthermore, for the barrier length N > 50 we can 
resort to a continuous description of the fork positions. The discrete master equation can then 
be rewritten as a bivariate Fokkcr-Planck equation [11,12] for the random walk of the two 
bubble-barrier boundaries characterized by the probability density function (PDF) P(x,y,t) 
of finding the boundaries at coordinates x — X/N and y = Y/N at a given rescaled time t 

with the dimcnsionless force / = N(u — 1)/(1 + u) and time rescaled by k(l + u)/2N 2 . 
Equation (J2J) is completed by the initial condition P(x, y, 0) = 6(x — xo)S(y — yo) and the 
reflecting boundary conditions (the bubbles in the soft zones are assumed to be open at all 
times) (d/dx — 2f)P(x, y, t)\ x=0 = (d/dy + 2f)P(x,y,t)\ y _ 1 = 0. Moreover, we impose the 
absorbing boundary condition P(x,x,t) = 0. This defines the vicious walker property [13], 
terminating the process when the two walkers meet. 

The dynamics of any number of vicious walkers in an arbitrary but common potential 
was studied previously [14], obtaining the solution by antisymmetrization of the single- walker 
evolutions (Green's functions). This method, however, already fails for two walkers in two 
different potentials. Of course, equation (JSJ) can be solved numerically as a (2+l)-dimcnsional 
partial differential equation, but it is clear that expressing the solution in terms of single- walker 
Green's functions would be very advantageous. It reduces the dimension of the problem by 1 
and even enables us to make exact analytical statements about some aspects of the solution. 

We show that for the particular problem at hand we can find the solution to equation 
by antisymmetrization of a suitable auxiliary single-walker Green's function. This is achieved 
by the similarity transform of the Fokker-Planck operator into a Hermitian operator (see 

reference [12], chapter 5) by using exp(=F/(z-a;o)) (|?r T 2 /^) cxp(±/(x - x )) = J^r - f 2 
which turns both Fokkcr-Planck operators into the same Hermitian one. Together with the 
freedom to choose the same boundary conditions for both walkers( 2 ) we end up with identical 
differential equations for both walkers and, thus, the antisymmetrization procedure can be 
used for the Hermitian version of the problem. This finally results in 

P(x, y, t) = e H*-y-*°+y°) (p(x, t\x )p(y, t\y ) - p(y, t\x )p(x, t\y )) , (3) 

where the auxiliary single- walker Green's function p(x, t\xo) satisfies the following ID equation 

= 0, (4) 

K=l/2=Fl/2 

( 2 )Each walker can reach only one boundary since the other one is blocked by the second walker. Therefore, 
we can choose the remaining boundary condition at will and, in particular, to be that of the second walker. 
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with the initial condition p(x, 0\xq) — 5{x — xq). Equation is solved by 

p(x, z\x ) = (e k]x - x °\ + K a e 2fc-fc|*-»o| + Ke Hx+x«) + ^^1] /2k [ K 2 e 2k - l] , (5) 

in the Laplace domain defined by g(z) = J °° g(t) e~>cp(—zt)dt. Here, k = \J ' z + f 2 and k = (k+ 
f)/(k— f) . To find the expression for the PDF P(x, y, t) we need the Laplace inverse of p(x, z) . 
Note that the latter can be expanded in eigenmodes through p(x, t) = Y^=o e *"Vn 0*0 

It is important to notice at this point that equation |J3J is not a physically meaningful ID 
Fokkcr-Planck equation (due to unphysical boundary conditions) and, indeed, it turns out that 
for / < one of the A's, determined by the transcendent equation (A + 2/ 2 ) sinh y/X + f 2 + 
2fy/X + f 2 cosh y/X + f 2 = 0, is actually positive, i.e., gives a solution exponentially growing 
in time. This type of behaviour is absent in the physical quantity P(x, y, t) due to the 
antisymmetrization procedure ©, yet the A > eigenmode plays an important role in the 
Kramers type asymptotics obtained for / <C — 1. Details are deferred to a future publication 
[10]. 

Results for the coalescence time and position. - The coalescence time PDF, 7r(f|xo, yo), 
can be found using the conservation of probability 

* /•! rv 

dt'ir{t')+ dy dxP(x, y, t) = 1, (6) 
Jo Jo 

where the two terms denote the probabilities that walkers have met before time t or have not, 
respectively. Then it follows from the Fokker-Planck equation J2J) after integration by parts, 

= / dxg(x,t), (7) 

y— x JO 

where g{x, t) is the joint PDF for coalescence time and meeting position of the two walkers. 
We are interested in nit) and its mean value (mean coalescence time) r = J Q tn(t)dt together 
with the PDF of the meeting position p(x) = / °° g(x,t)dt. By standard identities of the 
Laplace transforms we note that p{x) = g{x, z = + ) and r = — dir(z)/dz\ z=0+ . The related 
analysis becomes non-trivial as the auxiliary Green's function p(x, t\xo), as mentioned above, 
has a positive eigenvalue in the case / < 0. Details of the calculations using the theorem of 
residues and high precision numerics will be presented in reference [10]. 

Here, we discuss the final results for the coalescence time PDF ir(t) in figure|31 and for the 
mean coalescence time t and position PDF p(x) in figure 0] In all cases the initial condition 
used is xq = 0, yo = 1, i.e., the walkers start out at the reflecting walls or, in terms of the 
original model, the bubbles start from the boundaries of the soft zones. The coalescence 
time probability densities graphed in figure |21 for two different values of / (or u) compare the 
results of all three mentioned methods (Gillespie, master equation, and the semi-analytics) 
and show a sharp initial raise of rather peaky structures followed by an exponential decay 
for long times. The short-time increase stems from the initial separation of the two forks 
that first have to go through several random steps before having the chance to coalesce while 
the long time asymptotics is fixed by the boundary conditions. The agreement between the 
various computational approaches is very good. The locations of the maximum probability 
peaks correlate well with the termination points of sample trajectories depicted in figure|5]for 
the same parameters. 

The results for (i) p(x) and (ii) t calculated by the semi-analytical method are shown in 
figure 0] (i) The curves for the PDF p(x) of the coalescence position exhibit a pronounced 
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Fig. 3 - Comparison of the results for the bubble coalescence time PDF rr(t) obtained from the 
presented semi-analytical theory (eqs. by numerical inverse Laplace transform, the master 

equation, and the stochastic simulation (Gillespie), demonstrating excellent agreement. The values 
of the parameters are the same as in figure^ i.e., u — 0.98 (left) and u = 1.10 (right), and N = 50. 



crossover from a relatively sharply peaked form to an almost flat behaviour. The former 
occurs for large positive force /, corresponding to a strong drift toward a potential well, with 
negligible influence of the boundary conditions. In contrast, for large negative /, corresponding 
to a high barrier for coalescence, the insensitivity of p(x) to the position x can be explained in 
terms of a simple Arrhcnius argument: The probability of the walker to be at a position x is 
proportional to the Boltzmann weight, exp(— j3(f)(x)), where <f>(x) — — J F[x')dx' is the free 
energy corresponding to the force F(x) (see figure^. Then, the joint probability to have both 
walkers meet at the same position is given by the product exp(— 0[(f>L(x) + 4>r(x)]) ss const as 
the two walkers are in opposite linear potentials and the position dependence of the exponent 
cancels out. This simple picture necessarily breaks down close to the boundaries, (ii) The 
/-dependence of the mean first passage time r crosses over from the t ~ 1/ f behaviour typical 
for diffusion in a strong positive force pushing the two walkers together, to the exponential 
form t ~ exp(2|/|) of the associated Kramers problem. The former problem was studied in 
reference [15] by neglecting the boundaries and switching to the relative coordinate description 




Fig. 4 - Left: PDF p(x) of the coalescence position x for various values of the dimensionless force 
/. Right: Mean coalescence time r (in dimensionless, i.e., scaled, units) as a function of /. The two 
analytic asymptotic behaviours for large and small / are also shown. 
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which enables one to find the analytic result r = 1/(4/). For the Kramers problem (/ <C 
-1) the analytic solution for both p(x) = [1 - e~ 2 M x - e^K*-*)] \f\/(\f\ - 1) and r = 
e 2 '^'/[16/ 2 (|/| — 1)] can be found rather easily [10] by the expansion into the lowest two 
eigenmodes of p(x,t\xo). 

Conclusions. - In the present study we investigate two-bubble coalescence in a designed 
DNA construct consisting of two soft regions separated by a more stable barrier zone. We 
present a continuous semi-analytical theory yielding the coalescence time and position PDF's 
and show that the results agree well with the numerical results obtained with the discrete 
master equation and Gillespie stochastic simulation schemes. We note that for long barriers, 
whose sequence is arranged as random energy landscape, the propagation of the zipping forks 
may become subdiffusive in time [16] and the associated first passage time process is generated 
from the present results with the help of the subordination formalism [17]. The mathematical 
analysis reduces to the previously unaddressed case of two vicious walkers in opposite linear 
potentials. We demonstrate that its solution can be constructed by antisymmetrization of 
appropriate auxiliary one- variable densities. 
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